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1. INTRODUCTION 


1.1 Introduction 

The problems associated with fatigue were brought into the forefront of research by 
the explosive decompression and structural failure of the Aloha Airlines Flight 243 in 1988. 
The structural failure of this airplane has been attributed to debonding and multiple cracking 
along the longitudinal lap splice riveted joint in the fuselage. This crash created what may 
be termed as a minor “Structural Integrity Revolution” in the commercial transport 
industry. Major steps have been taken by the manufacturers, operators and authorities to 
improve the structural airworthiness of the aging fleet of airplanes. Notwithstanding this 
considerable effort there are still outstanding issues and concerns related to the formulation 
of Widespread Fatigue Damage which is believed to have been a contributing factor in the 
probable cause of the Aloha accident. The lesson from this accident was that Multiple-Site 
Damage (MSD) in “aging” aircraft can lead to extensive aircraft damage. A strong 
candidate in which MSD is highly probable to occur is the riveted lap joint. 

1.2 Background 

Riveted lap joints are used in an aircraft fuselage to join large skin sections. Among 
the many different types of joints, the single lap riveted joint is commonly used in aircraft 
construction. Joining introduces discontinuities (stress raisers) in the form of holes, 
changes in the load path due to lapping, and additional loads such as rivet bearing and 
bending moments. Because of these changes at the joint, local stresses are elevated in the 
structural component. Accurate estimations of these local stresses are needed to predict 
joint strength and fatigue life. 

Exhaustive studies on stress-concentration factors (SCF'sj for holes and notches in 
two-dimensional bodies subjected to a wide variety of loadings have been reported in the 
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literature [1,2]. Studies have also been made on three-dimensional stress-concentrations at 
circular holes in plates subjected to remote tension loads [3-6], A paper by Folias and 
Wang [6] provided a review of these previous solutions and presents a new series solution. 
The Folias and Wang solution covers a wide range of ratios of hole radius to plate 
thickness. The stress concentration at a hole in a plate subjected to bending was first 
presented by Neuber [4] using the Love-Kirchhoff thin plate theory [7], Reissner [8] 
rederived the plate solution including the effect of shear deformation and showed that 
Neuber’s solution was unconservative. Reissner’s SCF solution for bending loads is 
presented in terms of the Bessel function. Naghdi [9] extended Reissner’s analysis to 
elliptical holes using Mathieu’s functions. Rubayi and Sosropartono [10] conducted 3-D 
photoelastic measurements to verify Reissner’s circular hole and Naghdi’s elliptical hole 
solutions. Other analytical solutions are given in references [11, 12]. Information on the 
fatigue behavior of riveted joints has been derived mainly from investigations associated 
directly or indirectly with aircraft. Experimental tests are usually performed on single lap or 
butt type joints, made with aluminum alloy plate and rivets, and loaded in repeated tension 
[13-15], Results are reported in literature for remote loading, but very few papers 
considered [16, 17] 3-D effects for rivet loading in the hole. 

A wealth of data on stress concentration at cut-outs in plates subjected to remote 
tension, remote bending, or simulated pin loading have been reported in the literature. 
Shivakumar and Newman conducted exhaustive 3-D finite element analysis of plates with 
holes and developed 3-D stress concentration solutions. Results were reported in the form 
of simple equations and a computer program [16], The SCF along with the S-N diagram of 
the material may be adequate for designing riveted joints. This design is conservative, 
because a high factor of safety has to be used to account for the various unknowns. 

With the advent of powerful computers, it became possible to explore this field by- 
using the Unite element method to simulate real situations. The work of Ekvall [18] is one 



example. He developed a simple finite element model for the stress analysis of a joint to 
determine the local stress and strain at the fatigue critical location of a riveted lap or splice 
joint. Then, fatigue life predictions based on the local stress at the critical point were made 
using the effective stress-life law and the predicted fatigue and experimental lives were 
compared. In the simple finite element model he developed, the rivets were only modeled by 
three spring constants corresponding to the stiffness due to axial load, a shear load and 
bending moment applied to the rivet. The contact between the rivet and the plates was 
ignored. 

A numerical and experimental fatigue analysis of cold worked holes with clearance 
fit and interference fit fasteners was made by Rich and Impellizzeri [19]. the equivalent 
strain amplitudes were computed from an equation developed by the authors and then 
entered on a strain amplitude-life plot of constant-amplitude data from smooth specimens 
and a fatigue life determined. The total life predictions were in reasonable agreement with 
the test data. However, both the analytical model and the test specimen simulating a fighter 
aircraft wing skin, where it would attached to a spar cap, were designed as a single piece and 
the load transfer through fasteners was ignored. Swenson et al. [20] developed a finite 
element model to simulate crack growth in the spanwise lap joints of an aircraft wing, where 
the primary loading is parallel to the joint. In their model, each layer of the riveted joint was 
represented by a separate 2-D finite element mesh and the layers were connected by rivet 
elements which were modeled as springs. It was found that the predicted crack growth rate 
was much higher than the test data at the start of the crack growth, although the predicted 
crack growth rate approached the experimental result as 
the crack length approached half the joint overlap. 

One major unknown is the load transfer mechanism at the joint. The load transfer 
mechanism is much more complex than the simple supeiposition of various loadings that 
can be used. Further more such a problem has not been solved. The complexities are 



surface-to-surface contact between the rivet and the plates, the joint rotation due to non- 
axiality of loading and nonsymmetry of the configuration, and rivet clamp-up and 
interference. Because the problem is 3-D, the complexity is increased by one order of 
magnitude. Furthermore, the contact deformation is nonlinear, hence requires the solution 
of a variable BVP (boundary value problem). Analysis of riveted joints including these 
factors is important for the efficient design of joints, establishing the true factor of safety, 
and to verify the adequacy of the present design guidelines. 

1.3 Total Fatigue Life Prediction Models 

In the riveting process, a head on the protruding end of the shank is formed and the 
shank is deformed and expanded laterally to fill the rivet hole. In doing so it naturally grips 
the plate together. When such riveted joints are subjected to frequently repeated cyclical 
loads during service, the stress concentration effect may produce a fatigue crack and finally 
the riveted joints will fail from fatigue cracking even though the maximum applied stress is 
still within the range of the elastic behavior of the body. 

There are two philosophically different approaches to predict the total life of a 
component. The difference in these methods lies in the definition of a crack. In approach 
one, the total life is defined as the sum of fatigue crack initiation life and crack propagation 
life. The second approach is called the small crack theory[21, 22], This method assumes 
all materials have cracks either as a physical crack or as a defect. Defects causes cracks in 
the very first load cycle. So the crack is formed. Therefore, the total life is total number of 
load cycles required to make these micro-cracks to cause the fracture. The size of the 
micro-cracks is same as the material defect caused by material processing. Although the 
first approach has been used in the industry for a long time, the small crack theory is 
becoming attractive because the total could be predicted from computer models. The small 
crack theory requires a computed description ot the stress field at critical locations. The 



crack growth is calculated under those stress field. The micro-crack propagation is 
calculated under the influence of stress concentration. When the crack becomes one-tenth 
of a millimeter, the complete stress field will be used for crack propagation. Therefore, 
stress analysis of a joint including all joint complexities is critical to successful prediction of 
the total life of the joint using small crack theory. 

1.4 Rivet Clampup and Interference 

The riveting process consists of inserting the rivet in matching holes of the pieces to 
be joined and subsequently forming a head on the protruding end of the shank, the holes 
are generally 1/16 in. greater than the nominal diameter of undriven rivet. The head is 
formed by rapid forging with a pneumatic hammer or by continuous squeezing with a 
pressure riveter. The latter process is confined to use in shop practice, whereas pneumatic 
hammers are used in both shop and field riveting. In addition to forming the head, the 
diameter of the rivet is increased, resulting in a decreased hole clearance or the expansion of 
the hole (interference) [23], 

Most rivets are installed as hot rivets, but some shop rivets are driven cold. Both 
processes introduces clampup force and interference to the joint. 

During the riveting process the enclosed plies are drawn together with installation 
bolts and by the rivet equipment. As the rivet cools, it shrinks and squeezes the connected 
plies together. A residual clamping force or internal tension results in the rivet. The 
magnitude of the residual clamping force depends on the joint stiffness, critical installation 
conditions such as driving and finishing temperature, as well as driving pressure. 
Measurements have shown that hot driven rivets can develop clamping forces that approach 
the yield load of the rivet. Residual clamping forces are also observed in cold driven rivets. 
This results mainly from the elastic recovery of the gripped plies after the riveter, which 
squeezed the plies together during the riveting process, is removed. Generally, the clamping 



force in cold-formed rivets is small when compared with the clamping force in similar hot- 
driven rivets. The clamping force in the rivet is difficult to control, however a range of 
clampup force as a percentage of the rivet yield load can be assumed for analysis. 

The critical joint component in a lap joint subjected to repeated loading is not the 
fastener but the plate material. A severe decrease in the plate fatigue strength is apparent in 
unrestrained lap joints. The inherent bending deformations cause large stress ranges to 
occur at the discontinuities of the joint. The bending stress combines with the normal 
stress and results in high local stresses that reduce the fatigue strength of the lap joint. 

1.5 Problem Definition 

Two joint problems were considered; one is a classical pin joint (see Figure 1.1) to 
establish the basic solutions and the other is a two rivet single lap joint (see Figure 1 .2). All 
the dimensions are in millimeters. The lap joint is an idealized version of a riveted panel 
tested by Hartman [24], see Figure 1.3. The total life of this panel was measured [24] and it 
was predicted by the small crack theory by Newman et al [21,22], Newman used the results 
presented in this thesis for the prediction of the total life of the panel. 

The pin joint configuration is almost the same as the lap joint, but the out-of-plane 
load eccentricity is eliminated. Therefore, effect of clampup, interference, and friction 
should be same as the lap joint. 

The lap joint solutions should demonstrate the load eccentricity effects. A detailed 
finite element analysis of the joints are conducted and the results are presented in the thesis. 
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Figure 1.1. Pin joint. 



Figure 1.2. Two rivet single 
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B = 1 mm Rivet head (D) 


Figure 1.3. Experimental panel. 

1.6 Objectives of the Research 

Objectives of this research are as follows: 

1 . Conduct a detailed 3-D stress analysis of the pin joint and double row 
single lap rivet joint including nonlinear contact and large deformation. 

2. Identify the regions of contact and high stresses, and then establish stress 
concentration factors. 

3. Assess the effect of rivet clamp-up, rivet interference and friction on the 
local stress. 


1.7 Scope 

The report is presented in six sections. The first section presents introduction, 
background, definition of the problem and objectives of research. Section 2 presents the 
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description of finite element modeling of the two joints and modeling rivet clampup and 
interference effect. Also discussed in this section is the convergence criteria used for the 
non-linear analysis. Section 3 details the pin joint analysis for the neat fit, friction, clampup, 
interference and the combined case. The combined case is a combination of the rivet 
clampup, interference and friction. Section 4 covers the two rivet single lap joint analysis. 

In section 5 the neat fit case has been extended to elastic - plastic analysis to simulate a 
more realistic condition and the local stresses in the two rivet joint. Conclusions from the 


study are summarized in section 6.. 



2. FINITE ELEMENT MODELING OF RIVETED JOINT 


2.1 Introduction 

This section describes the finite element modeling of the rivet joint, contact, friction, 
clampup and interference. The general analysis procedure and convergence criteria 
are presented. 

2.2 Finite Element Analysis 

A stress analysis problem involves the differential equations of equilibrium and 
compatibility, together with the stress strain relationships and the boundary conditions. 
Analytical solutions to real life problems are seldom possible, and it is necessary, therefore, 
to employ a numerical method. 

A number of numerical stress analysis techniques are currently available, and their 
implementation is being greatly facilitated by the increasingly widespread availability of 
computers. The essential common feature of these methods is that the original problem, 
posed in terms of differential equations in the unknown continuous functions, is replaced by 
a formulation involving a set of algebraic equations in the discrete values of the unknowns at 
a finite number of points in the solid. In other words, the continuum model of the problem 
is approximated by a discrete model having a finite number of degrees of freedom. 

Of the numerical methods available the finite element method is the most widely 
used. The finite element method is a numerical procedure for obtaining solutions to many 
of the problems encountered in engineering analysis. It is impossible to document the exact 
origin of the finite element method because the basic concepts have evolved over a period of 
150 or more years. The method as we know it today is an outgrowth of several papers 
published in the 1950s that extended the matrix analysis of structures to continuum bodies. 
The space exploration of the 1 960s provided money for basic research, which placed the 
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method on a firm mathematical foundation and stimulated the development of multiple- 
purpose computer programs that implemented the method. The design of airplanes, 
missiles, space vehicles, and the like, provided application areas. Although the origin of the 
method is vague, its advantages are clear. The method is easily applied to irregular shaped 
objects composed of several different materials and having mixed boundary conditions. It 
is applicable to steady-state and time dependent problems as well as problems involving 
both geometric and material nonlinearity. 

The finite element method combines several mathematical concepts to produce a 
system of linear or nonlinear equations. The number of equations is usually very large, 
running to several thousand depending on the problem that is being solved, and requires the 
computational power of the computer. The method has little practical value if modem 
computers are not available. The basis of the method is the representation of a structure by 
an assemblage of subdivisions or finite elements as shown in Figure 2.1. These finite 
elements are considered to be connected at joints, called nodes or nodal points, at which the 
values of the unknowns (usually the displacements) are to be approximated. Successive 
finer discretization of the structure would lead to the exact solution. Therefore, it is likely 
that a moderately fine subdivision will provide a solution of acceptable accuracy. The 
computational effort required to obtain a solution will depend upon the number of degrees 
of freedom in the finite element model. In engineering practice a limit will be imposed on 
the degree of subdivision of the structure by the need to strike a balance between computing 
costs and solution accuracy. 



Nodal point 


Element 



Figure 2.1. Representation of a two-dimensional solid as an assemblage of 
triangular finite elements. 

Numerous commercial finite element analysis software packages are now available 
for simulating and solving complex engineering problems. One such code is ANSYS [25]. 
One of the advantage of ANSYS is its capability for geometric modeling and post- 
processing. Geometric modeling, analysis and results visualizations are all in the package. 
The analysis options include static, dynamic, material and geometric nonlinear analysis. In 
addition to having standard 1 -D, 2-D, 3-D elements, it has line to line and surface to surface 
contact elements. These elements are needed for the present analysis of riveted 
joint. 

2.3 Finite Element Modeling of Rivet Joint 

The finite element model of the rivet joint (refer to Figure 1.1 and 1.2) consists of 
three main components namely the top plate, bottom plate and the rivet. 

The plates and rivet are discretized using the SOLID45 3-D 8-Node Structural Solid 
element. The element is defined by 8 nodes having three degrees of freedom per node 
(translations in the nodal x, y, and z directions). The element may have any spatial 
orientation. The element has plasticity, stress stiffening, large deflection, and large strain 
capabilities. It can tolerate irregular shapes without much loss of accuracy. SOLID45 
element has compatible displacement shapes and are well --uited to model curved 
boundaries. 
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Contact occurs between the top plate and the bottom plate, the rivet and the plate 
holes, the rivet head and the plate. In the ANSYS program general contact is a boundary 
nonlinearity feature that permits surface-to-surface contact analysis with large deformations, 
contact and separation, coulomb friction sliding, and heat transfer. General contact is 
represented in the ANSYS program by following the position of points on one surface (the 
contact surface) relative to lines or areas of another surface (the target surface). The 
program uses contact elements to track the relative positions of the two surfaces. Contact 
elements are triangles, tetrahedron or pyramids, where the base is made up of nodes from 
the second surface (the target surface) and the remaining vertex is a node from the 1st 
surface, the contact surface. An analysis that incorporates general contact surfaces can 
easily require the use of hundreds or even thousands of contact elements. Fortunately, 
special features have been included in the ANSYS program to make generating and using 
these elements as efficient as possible. During solution, the program identifies those 
relatively few contact elements that are expected to affect the solution (i.e. those approaching 
contact or in contact). The remaining elements are temporarily ignored, producing null 
element stiffness matrices. As a result, an increase in the number of contact elements that 
are not in contact will not degrade run times as severely as would a similar increase 
involving other element types. The contact element used for the present problem is 
CONTAC49 3-D Point to Surface Contact. 

2.3.1 CONTAC49 Element Description 

CONTAC49 is a 5 node element that is intended for general contact analysis. In a 
general contact analysis, the area of contact between two or more bodies is generally not 
known in advance. In addition the finite element models of the contacting bodies are 
generated in such a way that precise node-to-node contact is neither achievable nor desirable 
when contact is established. The CONTAC49 element has the capability to represent 



general contact of models that are generated with arbitrary meshes. In other words, its use 
is not limited to known contact or node-to-node configurations. 

CONTAC49 is applicable to 3-D geometry. It may be applied to the contact of 
solid bodies or shells, to static or dynamic analyses, to problems with or without friction, 
and to flexible-to-flexibie or rigid-to-flexible body contact. 

Contact Kinematics 

Contact kinematics is concerned with the precise tracking of contact nodes and 
surfaces in order to define clear and unambiguous contact conditions. The primary aim is 
to delineate between open (i.e., not in contact) and closed (in contact) contact situations. 
This task is accomplished by various algorithms embedded in the CONTAC49 element. 

Contact and target definition 

With reference to the Figure 2.2, two potential contact surfaces are referred to as 
either the “target surface” or the “contact surface”. The target surface is represented by 
target nodes I, J, K and L, and the contact surface is represented by the contact node M. It 
is usually the case that many CONTAC49 elements will be needed to fully represent a 
realistic contact problem. 



Contact Surfaces and Nodes 



Figure 2.2. The CONTAC49 element configuration. 

Pinball Algorithm 

In simple terms, contact occurs whenever the contact node M penetrates the target 
surface (I, J, K, L). The first step in the determination of contact penetration is to make a 
distinction between near-field and far-field contact. Referring to 2-D transverse plane view 
for simplicity, (see Figure 2.3) shows several positions of a contact node with respect to the 
target surface. For CONTAC49 in 3-D the delimiting region is a sphere which is referred 
to as the “pinball”. When a contact node is outside the pinball an open contact condition is 
assumed, irrespective of whether or not the contact node is above or below the target, 
penetration can only occur once the contact node is inside the pinball. The radius of the 
pinball is internally fixed to be 50% greater than the two target surface diagonals. 

Pseudo Element Algorithm 

The next step in the determination of contact is to associate a single target to each 
contact node depending upon the position of the contact node in space. This is 
accomplished by establishing solid “pseudo elements" for each target surface as shown in 


Figure 2.4. A unique association is formed whenever contact node M is found within a 
target’s pseudo element. If a clear distinction is not made it is possible that contact 
“voids” or “overlaps” can appear. These voids and overlaps are unavoidable and are due 
to piecewise discretization of surfaces that are actually curved. These solid elements are 
temporarily formed each equilibrium iteration and provide a continuous mapping for each 
contact node that is in or nearly in contact with a target. The kinematic information that is 
needed to build these pseudo elements is stored in a global contact data base that is updated 
each equilibrium iteration. 
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Figure 2.3. Definition of Near-Field and Far-Field Contact. 
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Figure 2.4. Pseudo Element. 

Contact Gap And Projection 

The pinball and pseudo-element algorithms provide a one-to-one mapping between a 
contact node and a target. The final kinematic step is to determine the open gap or the gap 
penetration of the contact node on the target plane, along with the point of projection of the 
contact node. This is achieved by first modifying the target surface nodes to lie in a plane if 
they do not already, simplifying tangential surface calculations. In other words the warping 
of the target surface is ignored. In Figure 2.5 several coordinate systems are indicated. The 
global system is the usual X-Y-Z system. The next system is the natural s-t-n system of the 
planar target surface. Another coordinate system x-y-z is constructed from the natural s-t-n 
system in such a way that n and z directions are parallel. This enables straightforward 
tracking of the tangential contact motions. Finally a second Cartesian x e , y e , z c system is 
defined tor the element. Having defined the modified (unwarped) target surface and the 
various coordinate systems, the contact kinematics of gap and location are left to be defined. 
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With reference to Figure 2.5, the contact location (s*, t*) is computed by Newton’s 
iterative method based upon a normal projection of the contact node to the target plane. At 
the projected contact point a value of gap (g) is determined by the contact node’s location 
with respect to the target plane. Contact penetration is assumed to occur if the value of g is 
found to be negative, and the s* and t* projections are found to be in the natural space 
bounds of the target. For the later condition, the target surface is internally expanded based 
on the tolerance specified, thereby increasing the chances that a contact node will come into 
contact with the target plane. A positive gap value indicates an open contact condition. 



Contact Forces 

As explained above, contact is indicated when the contact node M penetrates the 
target surface defined by target nodes I, J, K, and L. This penetration is represented by the 
magnitude of the gap (g) and is a violation of compatibility. In order to satisfy contact 
compatibility, forces are developed in a direction normal (n-direction) to the target that will 
tend to reduce the penetration to an acceptable numerical level. In addition to compatibility 
forces, friction forces are developed in directions that are tangent to the target plane. The 
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normal and tangential friction forces that are described here are referenced to the local x-y-z 
system shown in Figure 2.5. 

Normal forces 

Two methods of satisfying contact compatibility are available for CONTAC49: a 
penalty method and combined penalty plus lagrange multiplier method. The penalty method 
approximately enforces compatibility by means of a contact stiffness (i.e., the penalty 
parameter). The combined approach satisfies compatibility to a user defined precision by 
the generation of additional contact forces (i.e., Lagrange forces). 

For the penalty method, 

if g < 0 
if g > 0 



where K n is the contact stiffness (real constant KN). 

For the combined method, the Lagrange multiplier component of force is computed 
locally (for each element) and iteratively. It is expressed as 
f n = min(0, K n g + A i+I ) 

Where : A j+ , = Lagrange multiplier force at iteration i + 1 

_ fX, +aK n g if |g| > e 

U if 1st < £ 

£ = user - defined compatibility tolerance (Input quantity TOLN 
on R command 

cr = an internally computed factor (a < 1) 

Friction forces 

The CONTAC49 element considers three friction models: frictionless, elastic 
coulomb friction, and rigid coulomb friction. The Coulomb friction representations require 
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the input of the coefficient of sliding friction (jT). Friction causes the tangential forces, as 
the contact nodes meets and moves along the target surface. 

For the frictionless case the tangential force is merely: 

f, = f, = o 



For elastic Coulomb friction it is necessary to calculate the tangential deformation of 
the contact node relative to the target. Figure 2.6 shows the total motion (u) of the contact 
node M along the target plane. It is seen that the total tangential displacement (T|) is 
represented by the projection of the total contact node motion to the unwarped plane of the 
target. Two projection points are mapped in the natural coordinates (s, t). The point (s*, t*) 
is the current projection position, and the tangential deformation is tracked from the point 
(s„*, t 0 *) that is associated with the previous converged solution (i.e., the previous time 
point). The deformation is first separated into \ and v components, such that 

n = (: n; + n -) - 




where: 


Tj x = component of r\ in the local x direction 
T] y = component of rj in the local y direction 


Next, the deformation is decomposed into elastic (or sticking) and sliding (or inelastic) 
components. 

n * = K + r \\ 

Vy = % + nl 

Related tangential forces are: 

/* = m 

fy = Krf 

where: K t = sticking stiffness 

It follows that the magnitude of the tangential forces is 

f, = (/„’+/,')* 

The stiffness and the load vector for the CONTAC49 element is given below 

jX} r =[0 0 q, 0 0 q 2 0 0 q 3 0 0 q 4 0 0 1] 

{^,} r =[q, 0 0 q 2 0 0 q 3 0 0 q 4 0 0 1 0 0] 

{V} 7 =[° q, 0 0 q 2 0 0 q 3 0 0 q 4 0 0 1 0] 

For the 4- node target, individual interpolates are 

^3 =--( 1+ -y‘)( 1 +^) 

<? 4 =- i -( 1 -. y *)( 1 + r *) 


In the noimal direction, the force applied to the contact node (M) is balanced by opposite 
forces applied to the target nodes; that is, 



f n.M - f n,I + f n,J + f n,K + f n,L = f n 

Similarly, in the tangential directions, 


f *,M - f x,I + f„J + f x .K + f x,L = f * 
f y,M = f y,I + f y,J + f y,K + f y,L = f y 


Using the interpolation vector above, the element load vector (i.e., the Newton-Raphson 
restoring forces ) is: 


pr} = /.R}+/.{ w .}+/>,} 

and the stiffness matrix is given by 


M= 


*. K }K f + ^ K [N x }{/v,y + {w,}{Ar,}’ 

[o] 


if sticking contact 

if sliding or frictionless contact 
if open contact 


2.4 Modeling of Clampup 

There are two ways of modeling the clampup of the rivet. One is to shrink the rivet 
by applying a differential temperature (cooling the rivet) to the rivet and giving it only an 
axial thermal coefficient of expansion/contraction. The other method is to make the rivet 
shank length (T rs ) smaller than the thickness (T p ) of the plate. The differential dimensions 
induces a clampup pressure by the rivet head (see Figure 2.7). Magnitude of clampup force 
depends on the value of (T p - T rs ). Larger the value of (T p - T rs ) larger the clampup force. A 
calibration study was conducted to obtain the relationship between the clampup force and 
the rivet shortening. For the current problem the second method was used because it offers 
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a better control over the axial force on the rivet. The practical range of clampup force varies 
between 10% to 35% of the rivet yield load. 

Clampup Simulation (T rs < T p ): Clampup force = 10% to 35% of yield load of the rivet 


Figure 2.7. Schematic of the clampup procedure. 

2.5 Modeling of Interference 

The rivet interference was modeled using a larger diameter rivet to fit the joint. 
Figure 2.8 describes the mechanics of the modeling. The rivet hole diameter is R h and the 
rivet diameter is R,.. When the rivet is inserted in the hole, the differential (R,. - R h ) 
introduces the interference stresses. If the value is negative, there is no interference. 
Magnitude of interference depends on the value of (R r - R h ). Practical values of interference 
ranged from 0.00635 to 0.019 mm. 



nW — it-' 


rs 




Interference Simulation (R h < R r ): Interference (R r -R h ) = 0.00635, 0.0127, 0.019 mm 



Figure 2.8. Schematic of the interference procedure. 

2.6 Modeling of The Combined Case 

The combined case was modeled by simultaneously making the rivet shank height 
less than the thickness ot the plate and by making the rivet shank diameter to be larger than 
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the plate hole diameter. 


2.7 Analysis Procedure 

The commercial finite element code ‘ANSYS’ was used. The displacement method 
of analysis was used. The linear solution was obtained by the frontal solver. Before 
solution, ANSYS automatically reorders the elements for a smaller wavefront (smaller the 
wavefront less the CPU time required for solution). The nonlinear solutions are obtained 
from the Newton-Raphson iterative algorithm. The analysis was conducted by incrementing 
the displacement and calculating the equilibrium condition and the associated stress-strain 
field. The analysis was continued till desired stress state or the loading was 
attained. 

2.8 Convergence Criteria 

The force convergence criteria was used to solve the problem. This is the most 
efficient convergence criteria for nonlinear finite element problems. Since both nonlinear 
geometry and changing status elements were used in the model the convergence criteria was 
slackened to avoid convergence difficulties. The convergence criteria was arrived at in an 
iterative manner, slackening the convergence criteria whenever convergence problems were 
encountered. 

The finite element discretization process yields a set of simultaneous equations: 

[*]{«} = {f“} 

where: 

[K] = coefficient matrix 

{a} = vector of unknown degree of freedom values 

{A'} = vector of applied loads 
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If the coefficient matrix is itself a function of the unknown DOF values (or their derivatives) 
then the above equation is nonlinear. The Newton-Raphson method is an iterative process 
of solving the nonlinear equations and can be written as: 


where: 

[AT/"] = Jacobian matrix (tangent matrix) 

i = subscript representing the current equilibrium iteration 
{/T / = vector of restoring loads corresponding to the element internal loads. 
{Z 7 "} — = residua] or out of balance load vector. 

In a structural analysis, [AT 7 ] is the tangent stiffness matrix, {«.} is the 

displacement vector and {F" r J is the restoring force vector calculated from element 
stresses. 

The iteration process described continues until convergence is achieved. 
Convergence is assumed when: 

|{*}|| < e s R„, 

where {/?} is the residual vector: 

{«}= 

liwih (S Rfy (Euclidean norm) 
e R = tolerance value 

Convergence, therefore, is obtained when size of the residual (disequilibrium) is less 
than a tolerance times a reference value. The default out of balance reference value 
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2.9 Summary 

With the advent of modem day computers and their ability to crunch numbers, finite 
element analysis has gained favor in the industry as an essential tool in their design process. 
ANSYS finite element code was used in this research project. The main reason being its 
capability to simulate contact between two bodies and its capability to do nonlinear analysis. 
Also a method of introducing rivet clampup and interference to the rivet joint was developed. 



3. PIN JOINT ANALYSIS 


3.1 Introduction 

In this section a classical pin joint was modeled using 3-D finite elements. This 
joint is loaded by remote tension and is restrained by a pin. The effects of clampup and 
interference on the stress distribution in the hole boundary is presented in this section. 

3.2 Joint Configuration 

Figure 3.1 shows the geometry i.e. configuration of the pin joint. The plate was 
square with the edge being 20 mm and thickness being 2 mm. The hole was located in the 
center of the plate with a radius R h of 1 .6 mm. The pin head had a diameter (d h = 3.2R h )of 
5. 1 mm and thickness of 1 .0 mm (0.6R h ). The radius to width and radius to edge distance 
is greater than 6, hence the joint configuration represents the infinite plate configuration. 
The global Cartesian coordinate system is represented by x, y, z. The pin is fixed at its 
center line and the plate is pulled by a uniform displacement u 0 in the x direction. The pin 
bearing load (P) is the integral of x- directional reaction at the edge x = 10 mm. The 
geometry and the loading are symmetric about y = 0 and z = 0 planes. Hence, only one- 
quarter of the joint (shown by the shaded region) was modeled by finite elements. 

3.3 Analysis Model 

The joint symmetry was exploited to reduce computational time. Figure 3.2 shows 
the one - quarter geometry of the joint. The plate was loaded at x = 10 mm with a uniform 
displacement u x and the axis of the pin was fixed in the x and y direction (that is u x = u v = 0 
for the axis of the pin). The plate and the pin was modeled using 8-noded brick elements, 
SOLID45 in the ANSYS code. The 3-D surface to surface contact elements were used to 
simulate contact between the pin and the plate hole and pinhead and the plate surface. The 



finite element model had 6912 SOLID45 elements (3456 elements each in plate and pin) 
and 1920 CONTAC49 elements. Figure 3.3 shows the finite element mesh. 

The other boundary conditions imposed on the model were u v = 0 on y = 0 plane 
and u z = 0 on z = 0 plane. These two boundary conditions simulate symmetric deformation 
of the joint. 

3.4 Analysis Cases 

There are two types of non-linearities that are expected in the model, viz., nonlinear 
contact boundary and large rotation. Therefore, large deformation and non-linear contact 
strategies are used in the analysis. A commercial code ANSYS 5.3 was used. The non- 
linearities were solved by a modified Newton-Raphson iteration algorithm. The Lagrange 
multiplier and penalty methods are used for contact modeling. The defined maximum 
gap/penetrations and contact stiffnesses are 0.0 IHs and 2000 N/mnf t about 3% of the 
elastic modulus of the plate material, which was within the recommended range) 
respectively. The parameter Hs is the smallest element size in the model, which was 1/6 
mm. The residual force convergence criteria was used at every node to establish the 
convergence of the non-linear solution. The relative error in the nodal residual forces was 
less than 0. 1% of total applied force as a convergence criteria. 

The analysis was conducted for four different cases that occur in the joint: neatfit, 
clampup, friction and interference, and a selected combined case. The neatfit represents the 
baseline solution. This case represents no surface - surface friction, no clampup and no 
interference. Analytical modeling of each of these parameters are explained in section 2 and 
is summarized in the following sections. The analysis was conducted by incrementally 
loading the joint to an applied remote stress of about 1 50 MPa or about U 0 = 0. 1 mm. 
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3.4.1 Elastic Friction 

Elastic friction values of 0.1, 0.2, 0.3 and 0.8 were chosen. The stress analysis was 
conducted for all four values. Interestingly, both location and magnitude of local stresses 
were unaltered. Hence, it was concluded that friction has no effect on local stresses. 

3.4.2 Pin Clamp-up 

As explained in the section 2.3 the pin clamp-up was introduced by changing the 
length of the pin shank. By shortening the pin’s shank length compared to the thickness of 
the plate introduces clamp-up force in the joint. A separate stress analysis was conducted to 
establish a relation between clamp-up force and pin shortening. This relationship was 
found to be linear (refer to Figure 3.4 ). The clampup equation was given by 
Clampup force, F c = 64,054 * AL 
where AL is the rivet shortening (T rs - T ) in mm. 

The amount of pin shortening for clamp-up force of 10%, 25%, 35% the pin yield 
force was calculated. These values were 7.64, 19.1 and 26.7 Jim respectively. The analysis 
was repeated for all three values of pin shortening. 

3.4.3 Pin Interference 

Pin interference was introduced by increasing the radius of the pin (Rj.) in relation to 
the hole radius (R h ). Three values of interference 2(R,.-R h ) were chosen, namely, 12.7, 25.4 
and 38. 1 |llm. These values bound the real values experienced in the aircraft industry. 

3.4.4 Combined Case 

The combined case is a combination of rivet clampup, interference and friction 
between contacting surfaces acting simultaneously. For the combined case a clampup of 
25% , an interference of 12.7 jam, and a coefficient of friction of 0.3 was used. 


37 



3.5 Results 


Results of the analysis conducted for various cases are represented in the following 
subsections. First, neat fit (zero surface-surface friction) results are presented. Then the 
effects of clamp-up and interference on local stresses are examined. The primary focus was 
on the maximum hoop stress on the hole boundary and the hoop stress at 90° to the z-axis. 
The 2nd case represents location of maximum hoop stress for open hole problems. All 
local stresses are normalized by the remote stress (aj as much as possible. The remote 
stress was calculated by dividing the total reaction at x = 10 mm edge by the area of cross- 
section ( 10x1 mm 2 ). 

3.5.1 Neat Fit Results 

3.5. 1.1 Deformed Shapes 

Figure 3.5 show the global and local to pin deformed shapes of the joint at a load 
level of 156 MPa. As can be seen, the pin looses contact with the plate from 0 = 0° to 90° 
and then it maintains contact and the hole is deformed into an elliptic shape. 

3.5.1.2 Contact Non-linearity 

The effect of contact non-linearity on local stresses was examined by analyzing the 
hoop stress at 0 = 90° on the hole boundary of the plates. Figure 3.6 shows the variation of 
g 9 with remote applied stress (a„ ) at various values of ‘z’ at 0 = 90°. As can be seen a 9 
varies linearly with G„. The same results are plotted as stress concentration factor (SCF = 

<7 0 l <TJ >n Figure 3.7. The SCF is maximum for the bottom plate at Z = 0 (about 7.36) and 
lowest at z/t =0.5 (6.2) at the top surface of the plate. 

Figure 3.8 shows a linear variation of membrane and bending stresses with o„. 
Membrane stress is the average stress through the thickness and bending stress is half the 
difference between the top and bottom surfaces of the bottom plate at 0 = 90°. As can be 
seen the bending component is negligible compared to the membrane stress. The 
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membrane SCF was 6.78 and bending SCF was 0.57. Therefore, for a smooth fit rivet joint, 
the local stress field varies linearly with the applied remote stress. 

3.5.1.3 Radial Stress Distribution at the Hole Boundary 

Figures 3.9 and 3.10 shows the contact (radial) stress distribution on the hole 
boundary for remote applied stress of 29 and 156 MPa respectively. Radial stress is 
maximum at 0 = 1 80° for the plate. The radial stress can be approximated by cosine 
function as shown in Figure 3. 1 1 (see thick solid curves). These cosine functions can be 
represented in the form of 

— = a cos" 0 

<T, 

Values of ‘a’ and ‘n’ for various levels are given by the following table. 


Table 3.1 Values of ' a ' and ’it' for various 
remote stresses 


MPa 

z/t 

a 

n 

29 

0.00 

- 9.40 

0.88 


0.25 

- 7.40 

0.87 

60 

0.00 

- 7.30 

0.77 


0.25 

- 7.20 

0.90 

92 

0.00 

- 6.80 

0.83 


0.25 

- 7.70 

1 .00 

124 

0.00 

- 6.60 

0.92 


0.25 

- 7.70 

1.12 

156 

0.00 

- 7.20 

1.02 


0.25 

- 8.10 

1.19 


The rivet contact angle is defined as the angle over which the radial compressive. This angle 
was found to be nearly 90° for -3t/8<z<3t/8. Note that contact stress is zero at 0 = 90° for 
most of the locations through the thickness. Results at z = t p /2 (comer location) may not be 
accurate because they are being affected by rivet head contacts. 
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3.5.1.4 Hoop Stress Distribution 

Figures 3. 12 to 3. 13 show the hoop stress distribution around the hole boundary for 
various remote loads. Hoop stress is compressive for 0 < 30° and tensile for 0 > 30°. Hoop 
stress is maximum at 0 = 90°, for all values of z. The major maximum o 9 occurred at z = 0 
(mid-plane). The maximum hoop stress occurs at 0 = 90° for all loads. Also the curves are 
almost invariant showing a linear relationship between the remote and local hoop stress. 

3.5.1.5 Hoop Stress Contour Plots 

Figure 3.14 shows the contours of hoop stress on z = 0 plane. As it was shown in 
Figures 3.12 and 3.13, the hoop stress is maximum at 0 = 90°. The maximum hoop stress 
location is valid all through thickness of the plate, except at the outer surfaces of the plate. 

3.5.2 Clampup Force 

Stress analysis was conducted for three values of clampup force, namely, 10%, 25% 
and 35% of the yield load of the rivet. Both hoop and radial stresses on the hole boundaries 
were examined. The important findings are presented. Maximum hoop stress occurred at 
the hole boundary at 0 = 90°. 

Figure 3.15 shows the variation of the hoop stress at 0 = 90° and z = 0 mm or. the 
hole boundary for various clampup forces. The Figure clearly demonstrates that the hoop 
stress at the hole boundary decreases with increase of clampup force. However the 
magnitude of decrease is not significant. It is of the order of 15, 30 and 45 MPa from the 
neatfit case for 10, 25, 35% clampup forces. Compared to the magnitude of local stresses, 
these values are not significant. Therefore, clampup effect can be ignored at high remote 
loads. 

Figure 3.16 shows the variation of normalized hoop stress at 0 = 90° and z = 0 for 
neatfit and three values of clampup forces. At low remote stresses, clampup has very 
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pronounced effect, but at higher loads, the clampup effect is not significant. A similar 
conclusion can be derived for combined clampup and friction (see Figure 3.17). 

3.5.3 Interference 

Figure 3.18 shows variation of a e at 0 = 90° and z = 0 with remote stress for three 
values of interference. The interference introduces high initial tensile stress on the hole 
boundary. Because of this the hoop stresses for interference cases is larger than neat fit 
case. With increase in remote stress, the local hoop stress increases at a slower rate than the 
neat fit case. Finally at higher loads rate of increase of hoop stress is same for both neat fit 
and interference cases. Therefore, interference results show two slopes with a transition 
region. The transition stress is unique for unique values of the interference. 

3.5.4 Combined Case 

Figure 3.19 shows a comparison of neat fit, 10% clampup, 25.4 mm interference 
and combined case hoop stress distribution. Combined case is 10% clampup, 25.4 pm 
interference and a surface to surface friction p = 0.3. The comparison show that the effect 
of friction and clampup is very small compared to the interference. Hence they can be 
neglected. The interference effect is large and it has to be properly accounted to predict 
fracture and/or life of a single rivet (pin) joint. 

3.6 Summary 

Conducted detailed 3-D stress analysis of single pin joint with friction, clampup. and 
interference. Practical ranges of clampup force (0 to 35 % of axial yield load of rivet) and 
interference ( 0 to 38. 1 pm) were used in the parametric study. Results concluded that 

( 1 ) Maximum stress was the hoop stress on the hole boundary and occurred at 
90" to the load axis. The hoop stress was maximum at mid-thickness of the 
plate. 
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(2) The contact angle was found to be nearly 180°. 

(3) Elastic friction had negligible effect on local stresses (hoop) and hence it can 
be ignored. 

(4) Clampup effect was dominant at low applied loads. Clampup decreases the 
local hoop stresses. But at high applied loads, clampup effect is small. 

(5) Interference was a major factor that impacted the local stresses (hoop stress) 
around the rivet hole. Interference introduces local tensile hoop stress at the 
rivet hole. This initial stress reduces the rate of increase of local stresses 
with remote loads. This causes the local hoop stresses to be lower than the 
neat fit results at high load levels. 

(6) Contact, friction and rivet clampup nonlinearities were confined to low axial 
loads. At high loads, the response is nearly linear. 
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Figure 3.1. Joint configuration. 
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Figure 3.2. One - quarter of the joint model. 
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Figure 3.3. Finite element model of the pin joint. 
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(a) Global deformation 
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Figure 3.5. Deformed shape at the pin-hole boundary. 
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z values 
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Figure 3.8. Variation of membrane and bending stress components with remote stress. 
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Figure 3.9. Radial stress distribution around the hole boundary (<j„ = 29 MPa). 
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Figure 3.10. Radial stress distribution around the hole boundary = 156 MPa). 
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0, Deg 

Figure 3.11. Cosine fit for radial stress distribution around the hole boundary (a„ = 29 MPa). 
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Figure 3.12. Hoop stress distribution around the hole boundary (a M = 29 MPa). 
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boundary (a M = 156 MPa). 
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Figure 3.14. Contour plot of the hoop stress in the plate (s oo “ 156.04 Mpa) 



Figure 3.15. Hoop stress vs remote stress at z = 0 and 0 = 90° for different clampup. 
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Figure 3.16. Stress concentration factors Vs remote stress for various ciampup forces. 
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Figure 3.17. Hoop stress vs remote stress at / = 0 and 0 = 90° for elastic friction and clampup. 



Figure 3.18. Hoop stress vs remote stress at /= 0 and 0 = 90° for various interference values. 




Figure 3.19. Comparison for the combined 


4. TWO RIVET SINGLE LAP JOINT 


4.1 Introduction 

In this section a single lap joint with two rows of rivets was modeled using 3-D 
finite elements. This joint is loaded by remote tension and is restrained by using anti- 
symmetric boundary conditions as explained in the following sections. The effects of 
clampup and interference on the stress distribution in the hole boundary is presented in this 
section. 


4.2 Joint Configuration 

The two-rivet plate joint configuration with all the geometrical parameters are shown 
in the Figure 4. 1 . This geometry represents a test configuration used in reference [31]. The 
plan and sectional view of a two rivet, single lap joint shown in Figure 4.2 was analyzed. 

The rivet shank was straight with 1.6 mm radius R r and 2 mm height. The rivet has two 
button heads. The head radius and depth were 2.55 mm and 1 mm respectively. The two 
plates, top and bottom were 145 mm long, 20 mm wide (w) and 1 mm thickness (t). The 
spacing (s) between the rivet centers was 20 mm. The edge distance was (E d ) 10 mm and 
the length L was 125 mm. 

The joint was loaded in tension with an uniform displacement U 0 . The average 
remote stress was G„. The two ends of the joint were supported laterally, to simulate the 

experimental condition, as shown in Figure 4.2. The Cartesian co-ordinate system selected 
in modeling is also shown in Figures 4. 1 and 4.2. The X, Y, Z represents the global co- 
ordinate system and x, y, z represents the local co-ordinate system used for plotting results. 





Figure 4.2. Joint configuration of double rivet single lap joint. 


4.3 Analysis Model 

The rivet joint problem was symmetric about y = 0 plane, hence only one half of the 
joint was modeled ( Figure 4.3). To further reduce the size of the problem, only one quarter 
of the joint was modeled. This involved imposition of complex constraint condition on X = 
0 plane. This constraint condition simulates the mirror image of deformation in the two 
quarters (left and right). The deformation boundary condition on x = 0 plane is a 
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combination of anti-symmetry and skew symmetry. The boundary conditions can be 
expressed by the following constraint equations: 


U(0, y, -z) = -U(0, y, z) 
V(0, y, -z) = V(0, y, z) 


W(0, y, -z) = -W(0, y, z) 


These boundary conditions reduce the finite element model to one quarter of the 
joint (see Figure 4.4). 



Figure 4.3. Sectional 3-D view showing cyclic anti-symmetry. 



Figure 4.4. One fourth of the model. 


Both the plates and the rivet are modeled with 8 noded solid elements (SOLID45 in 
ANSYS). Finer idealization was used at high stress gradient regions. Various views of the 
finite element idealization are shown in Figure 4.5. The hole region was divided into 
segments of 7.5° and the plate thickness was divided into four layers. The rivet shank and 
head are idealized to match the rivet hole and plate (although the mesh matching was not 
required for the analysis). Global, local, plan, and sectional views of rivet and plate models 
are shown in Figure 4.5. The model had 4776 elements and 6650 nodes. The possible 
contact surfaces in this model were: 

1) Between the rivet shank and top and bottom plate’s hole. 

2) Between the inner surfaces of rivet head and plates. 



3) Between the bottom surface of top plate and top surface of bottom plate. 

These contact surfaces were modeled using 5-noded 3-D surface to surface contact 
elements represented by CONTAC49 in ANSYS code. Each target node has the possibility 
of contacting four elements. The model contained 2588 contact elements. 



Top & bottom plates 
and contact surface 


Figure 4.5. Various views of the rivet joint finite element model. 

The loading imposed on the model was a uniform displacement ‘U 0 ’ at X = L. In 
summary, the following boundary conditions are applied on the model: 

(1) Symmetry on y=0 plane. 

(2) Constraint equations on nodes at x=0 plane. 

(3) Uy = 0 at x = L, y = 0 and z = 0 (for restricting rigid body motion). 

(4) The loading (displacement) u = u 0 was imposed at x = L ( 1 25 mm) plane. 


o/ 


Aluminum alloy 2024-T3 Alclad properties were used in the analysis. The rivets are 
2024 (type DD) aluminum. The elastic modulus E = 68,950 MPa, Poisson’s ratio v = 

0.3, yield strength of 270 MPa and ultimate strength of 270 MPa. 

4.4 Analysis Cases 

Two types of non-linearities were expected in the model, viz., nonlinear contact 
boundary and large rotation. Therefore, large deformation and non-linear contact strategy 
were used in the analysis. A commercial code ANSYS 5.3 was used. The nonlinearities 
were modeled by modified Newton-Raphson iteration algorithm. The Lagrange multiplier 
and penalty method were used for contact modeling. The defined gap/penetrations and 
contact stiffnesses are about 0.0 lHs and 2000 N/ram 2 (about 3% of the elastic modulus of 
the plate material, which was within the recommended range) respectively. But for 
interference cases gap/penetration value used was 0.025Hs. Where Hs was the smallest 
element size in the model, which was 0.25 mm. The residual force convergence criteria was 
used at every node to establish the convergence of the non-linear solution. Relative error in 
the nodal residual forces was less than 0. 1% to 1% of total applied force as a convergence 
criteria . 

The analysis was conducted for three different complexities that occur in the joint. 
They are friction between contacting surfaces, rivet clamp-up, and rivet interference. 
Analytical modeling of each of these parameters is explained in section 2 and is summarized 
in the following sections. The analysis was conducted by incrementally loading the joint to 
an applied remote load of about 130 MPa or about U 0 = 0.3 mm. 
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4.4.1 Friction 


Friction between the contact surface was modeled as elastic coulomb friction. The 
surface tangent stiffness KT was selected to be KN/100, where KN was the normal contact 
stiffness. The tangential friction force at the contacted nodes was the product of friction 
coefficients and the normal force. The sliding friction coefficients used in the analysis were 
0, 0.3 and 0.8. The friction coefficient value of zero represents the smooth contact. 

4.4.2 Rivet Clamp-up 

As explained in section 2.3 the rivet clamp-up was introduced by changing the 
length of the rivet shank. By shortening the rivet length compared to the thickness of the 
two plates clamp-up force was introduced. A separate stress analysis was conducted to 
establish a relation between clamp-up force and rivet shortening. This relationship was 
found to be linear (refer to Figure 3.4). The clampup equation was given by 
Clampup force, F c = 64,054 * AL 
where AL is the rivet shortening (T rs - T ) in mm 

The amount of rivet shortening for clamp-up force of 10%, 25%, 35% rivet yield 
force was calculated. These values were 7.64, 19.1 and 26.7 pm respectively. The analysis 
was conducted for all these values of rivet shortening. 

4.4.3 Rivet Interference 

Rivet interference was introduced by increasing the radius of the rivet (R T ) in relation 
to the hole radius (R h ). Three values of interference 2(R,.-R h ) chosen were 12.7, 

25.4 and 38. 1 pm. These values bound the real values experienced in the aircraft industry. 


4.5 Results 

Results of the analysis conducted for various cases are represented in this section. 
First, neat fit (zero surface-surface friction) results are presented. Then the effects of 
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friction, clamp-up and interference on local stresses were examined. The primary focus was 
on the maximum hoop stress on the hole boundary and the hoop stress at 90° to the x-axis. 
The second case is where the hoop stress is maximum for open hole problems. 

4.5.1 Neat Fit Results 

4.5.1. 1 Deformed Shapes 

Figure 4.7 shows the sectional view of the deformed shape of the model. Notice 
that the two plates slided one past the other at the left end of the model. This is the true 
deformation expected if the complete half model was analyzed. This deformation pattern 
confirms the approximation of the boundary conditions imposed on X = 0 plane. The 
close-up view of the model at the rivet is shown in Figure 4.8. Notice that the rivet is in 
contact with top plate on the right side and bottom plate left side. The right side of bottom 
plate and the left side of top plate are separated from the rivet shank. The two plates are 
separated from each other on the loaded side. In the 3-D model like this problem, all 
separation and contact are only partial and they change along the polar (0) direction of the 
rivet. Figure 4.9 shows the superposed view of the two quarters of the model. This depicts 
the local joint rotation and deformation. 

4.5.1.2 Contact Nonlinearity 

The effect of contact nonlinearity on local stresses was examined by analyzing the 
hoop stress at 0 = 90° on the hole boundary of the plates. Figure 4. 10 shows the variation 
of G 0 with remote applied stress (a„ ) at various values of ‘z’ at 0 = 90°. Except at low 
values of applied stress, < 30 MPa, o e varies linearly with g„. The same results are 
plotted as stress concentration factor (SCF = o e / oj in Figure 4.11. Again, the SCF is 
maximum for the bottom plate at Z = 0, which is about 5.7. 

Figure 4. 12 shows almost linear variation of membrane and bending stresses with 
Membrane stress is the average stress through the thickness and bending stress is half 
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the difference between the top and bottom surfaces of the bottom plate at 9 = 90°. Both 
membrane and bending stresses vary linearly through the thickness. The membrane SCFs 
were 4.75 and 3.4 for the bottom and top plates whereas the bending SCF were 1 .2 and 0.6 
for the bottom and top plates respectively. 

4.5.1.3 Radial Stress Distribution at the Hole Boundary 

Figures 4. 13 and 4. 14 shows the contact (radial) stress distribution on the hole 

boundary for remote applied stress of 45 and 84 MPa respectively. Radial stress is 

maximum at 0 = 1 80° for the bottom plate and 0 = 0° for the top plate (not shown). The 

radial stress can be approximated by cosine function (see thick solid curves). These cosine 

functions are of the form shown below; 

ov „ 

— - = a cos 6 

o- 

Note that contact stress is zero at 9 = 90° for most of the locations through the thickness. 
Results at z = -1.0 (corner location) may not be accurate because they are being affected by 
rivet head contacts. 

4.5.1.4 Hoop Stress Distribution at the Hole Boundary 

Figures 4.15 and 4.16 shows hoop stress distribution on the hole boundary for the 
applied stress of 45 and 84 MPa respectively. The hoop stress peaks at 0 = 90° for values 
of z through the thickness. The maximum stress is at z = 0.0, which is about 5.7 times that 
of the remote stress. 

4.5.1.5 Hoop Stress Contour Plots 

Figure 4.17 shows the hoop stress contour plots on both top and bottom plates. 

The plate is sectioned at x = 0 plane to examine the location of maximum sq through the 
width. Results on Figure 4. 1 7 concludes that for a smooth fit case, the location of the 
maximum hoop stress is on the hole boundary and at z = 0.0. 
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In all the results presented so far, hoop stress on the hole boundary is highest and it 
occurs in the bottom plate at the interface between the two plates. 

4.5.2 Elastic Friction 

Stress analysis was conducted for elastic friction values of 0.3 and 0.8 for applied 
remote stress varying from 0 to about 140 MPa. Both radial and hoop stress distribution 
along the hole surface and through the thickness were examined for all cases and compared 
with neatfit results. Comparison of results showed very little effect of friction on magnitude 
and location of peak values. Tensile stresses in the bottom plate were higher than that in the 
top plate. Therefore, only results for the bottom plate are presented for this case and all 
other remaining cases. Figure 4.18 and 4. 19 shows hoop stress variation along the hole 
surface at various plate thickness locations for p = 0.3 and p = 0.8 respectively. For both 
cases a e is maximum at 9 = 90° and z = 0. 

Figure 4.20 shows variation of c e at 9 = 90° with remote stress ( ) for p = 0, 0.3 

and 0.8. Results of p =0 and p = 0.3 were almost identical and at low values of c^,, p = 

0.8 stresses were smaller than p = 0 case, but for > 55 MPa, all three curves merged. 
Therefore, elastic friction alone has very little effect or no effect on joint stresses. 


4.5.3 Clampup Force 

Stress analysis was conducted for three values of clampup force, namely, 10%, 25% 
and 35% of the yield load of the rivet. Both hoop and radial stresses on the hole boundaries 
were examined. The important findings are presented. Maximum hoop stress occurred at 
fhe hole boundary except tor very’ small value of remote applied stress. At very' low' applied 
stress, the maximum hoop stress occurred in the interior of the piate, this may be because of 
pinching effect caused by clamping. 
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Figures 4.21 to 4.23 shows the a 9 stress distribution on the hole boundary for 10% 
clampup force at remote stresses of 45.7, 94.5 and 143.7 MPa respectively. Maximum 
hoop stress location changed from 75° at low load levels to 90° at high loads. Again the two 
plate’s interface location was highly stressed region in the plates. Table 4. 1 summarizes the 
magnitude and location of maximum hoop stress for all three clampup forces and various 
applied stress values. The location of maximum stress is at 75° for both 25% and 35% 
clampup force for higher values of applied stress. 

The table also lists the hoop stress values at z = 0 and 0 = 90°, and the difference 
compared to the maximum stress. The difference increases with increased clampup force 
but decreases with increase of applied stress. For 10% clampup force, one can neglect the 
difference for > 50MPa. Knowing cr 9 at 0 = 90° and z = 0, for any condition, one can 
calculate back the maximum value of the hoop stress and its location using the table 4. 1 
results. 

Figure 4.24 shows variation of o 9 at 0 = 90° and z = 0 for the three values of 
clampup forces and for various a „ loads. The Figure clearly demonstrates that, as expected 
the hoop stress at the hole boundary also decreases with increase of clampup force. The 
relationship between the clampup force and the hoop stress can be expressed by a simple 
equation: 

Clampup — (^8 ^Neatfit ^lC^ 

where C 0 is the clampup stress, defined as clampup force divided by clamping area of the 
plate. In the equation, all stresses are expressed in N/mm 2 (MPa). At low applied stresses, 
clampup has nonlinear effect on the local stresses. 

Figures 4.25 and 4.26 show the effects ot clampun on membrane and bending 
components of the stress on the hole boundary at 0 = 90°. Membrane stresses decreases 
with increase in clampup, where as the bending stresses increases. 
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Figure 4.27 shows the variation of with radial distance from the hole 
boundary. The clampup effect is confined to clamped region, beyond that the hoop stress 
are unaffected by the clampup; particularly for higher remote stress (cr OT = 140 MPa). 
However at lower load levels, say 65 MPa, the influence of clampup goes slightly beyond 
the clampup region. The magnitude of change is minimal. Hence it terms of crack initiation 
and crack growth, clampup can increase the crack initiation life but once the crack has 
grown beyond the rivet head, clampup has no effect on crack growth. 

4.5.4 Clampup and Elastic Friction 

This analysis was conducted with the intention that the clampup force would increase the 
normal force between the plate surfaces and hence the friction effect (the tangential friction 
force is equal to the product of the normal force and the coefficient friction). Figure 4.28 
shows that friction still does not have any effect on the hoop stress distribution on the hole 
boundary. Hence we can conclude that friction has no effect on stress distribution in the 
joint. 

4.5.5 Interference 

Figure 4.29 and 4.30 show the hoop stress distribution around the hole boundary for an 
interference value of 12.7 mm and remote stress value of 48.7 and 78.7 MPa respectively. 

As can be seen from Figure 4.29 the maximum hoop stress occurs at 0 = 67.5° and as the 
remote load increases the location of maximum hoop stress shifts to 0 = 75° at remote load 
of 78.7 MPa (Figure 4.30). This trend can be also be inferred from table 4.2 which lists 
location of maximum hoop stress for all interference values. The maximum hoop stress 
location varies from 135° (for low load levels) to about 75° for high load levels. 

Figure 4.31 and 4.32 show the radial stress distribution around the hole boundary. 
As can be observed the radial stress became negative before 90° for all values of z. This is 
because of the interference effect. Since the rivet shank diameter is greater than that of the 
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plate hole, the rivet is in contact with more than quarter of the hole boundary and this 
induces compressive radial stresses on the plate hole even for 0 < 90°. 

The effect of interference on hoop stress at z = 0 and 0 = 90° is shown in Figure 
4.33 and the results are compared with neatfit case. Notice there are three distinct regions 
behaving differently. At low applied stresses, the local hoop stresses (g 9 ) increases with 
increase in interference value and it is larger than the neat fit case. At very large applied 
stress, local ct 9 for the neat fit case is larger than all the interference cases. In between these 
two regions, the results cross each other. However, interference introduces initial tensile 
stress on the hole boundary and the rate of increase of a 9 decreases with the amount 
interference. The slope of g 9 curve decreases from neatfit to 12.7, 25.4, and 38. 1 pm. All 
four curves cross each other at different applied stress values. The trend of the curves 
indicate that eventually ( at large g„ ) they will all be parallel to each other with neatfit case 
having highest g 9 and the largest interference having the lowest o 9 . From these results we 
can conclude: 

1 . Because interference case introduces initial tensile stresses at = 0, the 
alternating component ot the local stress due to alternating <?„, is lower than the neat fit 
results. 

2. Slope of the G 0 vs is smaller with increased interference. This in effect 
reduces the alternating component of the local stress c 9 . 

4.6 Summary 

Conducted detailed 3-D stress analysis of single lap joint with two rows of rivet 
including the effect of friction, clampup, and interference. Practical ranges of clampup force 
(0 to 35 % ot axial yield load of rivet] and interference ( 0 to 38. 1 pm) were used in the 
parametric study. Results concluded that: 



( I ) For the neat fit case hoop stress on the hole boundary was the highest and it 
occurred in the bottom plate at the interface between the top and bottom 
plates at 90° to the load axis. The maximum stress concentration factor was 
found to be 5.7. Contact nonlinearity was confined to low levels of applied 
load (0„.< 30 MPa). At higher applied loads the response is almost linear. 
The radial stress distribution on the hole boundary can be represented by a 

G r n 

cosine function of type = Cl COS 9 . Values of ‘a’ and ‘n’ depend 

< 7*0 

on the through the thickness location on the hole boundary. 

(3) Elastic friction has negligible effect on local stresses (hoop) and hence it can 
be ignored. 

(4) For the clampup case, again the maximum stresses occurred on the hole 
boundary and in the bottom plate. The clampup force did decrease the 
stresses at the critical location and its effect can be summarized in the form 
of a simple linear equation. For 10% clampup level, maximum hoop 
stress location changed from 75° at low load levels to 90° for high load 
levels. But for 25% and 35% clampup it remains at 75° for almost all load 
levels. The Clampup effect was confined to the region under the rivet head. 
Beyond the rivet head it’s effects were significantly reduced. 

(5) Interference is major factor that impacts the local stresses (hoop stress) 
around the rivet hole. Interference reduces the rate of increase of local 
stresses with remote loads. This causes the local (hoop) stresses to be 
lower than the neat fit results at high load levels. 
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Table 4.1. Maximum hoop stress location in bottom plate for different clampup forces 
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117.6 572.6 1.6 75.0 0.0 535.9 

143.4 705.6 L6 82.5 0.0 676.4 



Table 4.2. Maximum hoop stress location in botto m plate for different interference 

Remote Maximum Coordinates Stress %difference 

Stress Hoop stress r g z at z=0,q=90 

Interference = 12.7 mm 
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Figure 3.9. Deformed shape of the joint (super imposed image). 
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Figure 4. 10. Hoop stress vs remote stress for various z values at 0 = 90°. 
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Figure 4.12. Mean and bending stress for bottom plate. 
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Figure 4.13. Radial stress distribution around the hole boundary (s°o = 45.27 Mpa). 
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Figure 4.14. Radial stress distribution around the hole boundary (s°o = 84.26 Mpa). 
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Figure 4.15. Hoop stress distribution around the hole boundary (o M = 45.27 MPa). 
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Figure 4.17. Hoop stress contour for a section at 0 = 90° 




Figure 4.18. Hoop stress distribution at hole boundary (p = 0.3). 



Figure 4.19. Hoop stress distribution at hole boundary (|i = 0.8). 
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Figure 4.20. Hoop stress vs remote stress for various values of p. 




the hole boundary (cr^ = 45.7) 




Figure 4.22. Hoop stress distribution around the hole boundary (a^ = 94.5). 
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Figure 4.23. Hoop stress distribution around the hole boundary = 143.7). 
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Figure 4.24. Hoop stress vs remote stress at z = 0 and 0 = 90° for different clampup 



a a cl 

3 3 3 



o 

o 

© 

o 

o 

o 

o 

o 

o 

© 


o 

CO 

<o 

<3 * 

CM 


CM 

■ 


b 


98 


Figure 4.25. Membrane stress for bottom plate at 0 = 90°. 
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Figure 4.26. Bending stress for bottom plate at 0 = 90' 
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Figure 4.27(b). Variation of hoop stress in radial direction at 0 = 90° and z 
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Figure 4.28. Clampup + friction @ 0 = 90 and z = 0 of bottom plate. 
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Figure 4.29. Hoop stress distribution around the hole boundary (o M = 48.7 MPa). 
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Figure 4.30. Hoop stress distribution around the hole boundary (cr,*, = 78.7 MPa). 
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Figure 4.31. Radial stress distribution around the hole boundary (ct^ = 48.7 MPa) 
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Figure 4.32. Radial stress distribution around the hole boundary (o^ = 78.7 MPa). 
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Figure 4.34. Membrane stress vs remote stress at 0 = 90°. 
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igure 4.35. Bending stress vs remote stress at 0 = 90°. 


5. ELASTIC - PLASTIC ANALYSIS OF TWO - RIVET SINGLE LAP JOINT 


5.1 Introduction 

In the previous section, elastic analysis of the two - rivet joint was conducted. The 
local stresses at the rivet hole have exceeded the material yield stress. To investigate the 
redistribution of stresses due to plastic deformation, an elastic plastic analysis of the neat fit 
case was conducted. In this section results of this analysis are presented. The material was 
modeled as multi-linear stress - strain response, with an yield strength of 270 MPa and 
ultimate strength of 425 MPa. The joint configuration shown in Figure 1.2 and finite 
element mesh used is shown in Figure 4.5. 

5.2 Material modeling 

In this analysis the rivet material was assumed to be elastic and the plate material to 
be elastic - plastic, because the plastic deformation of the plate is of interest. The multilinear 
isotropic hardening stress - strain response was used to model the plate material. The stress 
- strain curve is shown in Figure 5.1. The von Mises yield criteria and incremental plasticity 
theory based on the associated flow rule was used in the analysis. The von Mises criteria 
for a three dimensional stress state in terms of the six Cartesian stress components is given 
by 




Yielding would occur when, 



where, <7 eff is the effective (von Mises) stress and cr ts is the current uniaxial yield stress of 
the material. Initially, a ys is the material yield stress (270 MPa) and as the material 
plastically deforms, the yield stress changes for isotropic hardening. In isotropic hardening, 
the yield surface expands in size as the plastic strains develop as shown in 
Figure 5.2. 

5.3 Analysis 

The plate was loaded incrementally by remote axial displacement u 0 (see Figure 
4.2). To follow the effective stress - strain response as closely as possible, small load 
increments were chosen. The initial displacement increment was 0.01 mm (3.92 MPa) and 
the subsequent displacement increments varied between 0.01 and 0.02 mm. The analysis 
was carried to a maximum remote displacement of about 0.22 mm. The force convergence 
criteria was used. That is the maximum out of balance residual force is about 1 to 10% the 
total incremental load for that load step. 

The analysis was conducted using a commercial finite element code ANSYS. The 
full Newton-Raphson method with updated tangent stiffness matrix and adaptive descent 
was used to solve this non-linear problem. Each analysis case took 9 CPU hours on DEC 
Alpha 255 computer. 

5.4 Results and Discussions 

The analysis was conducted for only one case, that is the neat fit ( no friction, no 
clampup and no interference) joint subjected to remote tensile loading. Only important 
results are presented here. 

Figure 5.3 show the progressive development of the plastic zone with remote 
loading on the hole boundary of the bottom plate. As can be seen yielding first occurs at 9 
= 1 80", where the rivet bears against the hole and then it progresses towards 0 = 90". The 



plate material starts to yield at around = 45.5 MPa (ux = 0. 1) at z = 0 and 0 = 1 80°. 

This yielding is due to the high compressive bearing stresses developed between the rivet 
and the plate. However compressive stresses are not critical in causing crack initiation or 
propagation. By the time remote loading reached 55.1 MPa, nearly half the rivet hole had 
yielded and at = 97.7 MPa, nearly 75% of the hole had yielded. At 0 = 90°, the yielding 
was due to high tensile stresses. Hence, the crack initiation and propagation is highly 
probable at 0 = 90°. 

Figures 5.4 through 5.6 show the variation of hoop stress around the hole boundary 
for three values of a„.. All stresses were normalized by the remote stress (a^.). With 
increasing load levels the stress concentration drops, for example from 5.4 (elastic) to 
around 3.45 (plastic) for <J„. = 97.7 MPa. Also, the peak tensile stresses are more evenly 
spreadout. Table 5. 1 lists the maximum hoop stress and its location and comparison of 
maximum hoop stress with a 9 at 0 = 90° and z = 0. In the elastic region the maximum 
stress location was at 0 = 90° and z = 0. After the material yields the maximum hoop stress 
location shifts to 75°. The difference between maximum hoop stress and a e at 0 = 90° and 
z = 0 is less than 1 .5 % except at g m . = 50.2 MPa. Therefore, we can conclude that ct 9 at 0 
= 90° and z = 0 are same the maximum hoop stress. 

Figure 5.7 shows the plot of o 9 at 0 = 90° and z = 0 versus the remote stress (c«, ) 
for elastic and elastic-plastic cases. The two results are almost same till = 45 MPa, and 
then they start diverging. The plastic analysis yields lower values of a e than the elastic 
analysis and a 9 reaches a plateau at about 320 MPa (to almost a + (a ult - cr )/3), The 
variation of membrane and bending components of hoop stress at 0 = 90° and z = 0 are 
shown in Figure 5.8. Both elastic and elastic-plastic solutions are presented. As stated 
previously, the elastic and elastic plastic analysis results starts diverging at = 45 MPa. 



The plastic membrane stress continue to increase and reaches a plateau. Whereas the stress 
decreases with and finally reaches almost a zero value. Therefore, we can conclude that 
the local bending stresses can be neglected at large remote stress. 

5.5 Summary 

Elastic-plastic analysis of neat fit two-rivet lap joint was conducted. The material 
was assumed to be multi-linear stress-strain response, Von Mises yielding and isotropic 
hardening. The incremental elastic-plastic analysis was used. The analysis was carried out 
to a remote stress of about 100 MPa. Hoop stresses at the critical hole at the joint was 
examined. The results showed that plastic deformation reduces stress concentration factor 
from 5.4 to 3.45, the hoop stress at 9 = 90° and z = 0 is nearly same as the maximum hoop 
stress and the bending stress reduce to nearly zero at large plastic deformation at the joint. 
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- strain response 
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Figure 5.2. Initial and subsequent yield surface for isotropic hardening 
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Figure 5.3. Progressive plastic - zone with loading in the bottom plate of the joint 



Elastic - Plastic anal 



Figure 5.4. Variation of hoop stress along 0 in the bottom plate at o, = 50.2 MPa 
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Figure 5.5. Variation of hoop stress along 0 in the bottom plate at = 74.3 MPa 
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Figure 5.6. Variation of hoop stress along 0 in the bottom plate at a„= 97.7 MPa 



Figure 5.7. Variation of hoop stress at 0 = 90° with cr, 
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Figure 5.8. Variation of membrane and bending component stresses with remote stress 
(at 0= 90° and z = 0) in the bottom plate 



6. CONCLUDING REMARKS 


Conducted a detailed stress analysis of pin joint and two-rivet single lap joint. 

These joint configurations are representative of lap-splice joints used in aircraft. The lap 
joint includes the load eccentricity effect where as the pin joint does not include this effect. 
The analysis included the effect of plate-rivet surface-to-surface friction, rivet clampup, and 
rivet interference effects. The analysis was conducted using 3-D, 8-noded isoparametric 
finite elements. The contact was modeled by 5-noded surface-to-surface contact elements. 
This element had capability of developing contact and separation between predetermined 
surfaces. Simple rivet contraction and expansion models were developed and used to 
simulate rivet clampup and interference. The range of clampup force used was 0 to 35% 
rivet yield load; rivet interference used was 12.5 to 38. 1 pm and the friction varied from 0.0 
to 0.8. A commercial finite element code was used for the analysis of the joints. Elastic 
analysis was conducted on the pin joint and elastic and elastic -plastic analysis were 
conducted on two-rivet single lap joint. The primary focus of results analysis was on hoop 
and radial contact stresses at the hole boundary. Results of this study yielded following 
conclusions. 

Elastic Analysis of Two-Rivet Joint 
Pin Joint 

( 1 ) Maximum stress was the hoop stress on the hole boundary and occurred at 
90° to the load axis. The hoop stress was maximum at mid-thickness of the 
plate. 

(2) The contact angle was found to be nearly 180°. 
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( 3 ) 


Elastic friction had negligible effect on local stresses (hoop) and hence it can 
be ignored. 

(4) Clampup effect was dominant at low applied loads. Clampup decreases the 
local hoop stresses. But at high applied loads, clampup effect is small. 

(5) Interference was a major factor that impacted the local stresses (hoop stress) 
around the rivet hole. Interference introduces local tensile hoop stress at the 
rivet hole. This initial stress reduces the rate of increase of local stresses 
with remote loads. This causes the local hoop stresses to be lower than the 
neat fit results at high load levels. 

(6) Contact, friction and rivet clampup nonlinearities were confined to low axial 
loads. At high loads, the response is nearly linear. 

Two-Rivet Lap Joint 

(1) For the neat fit case hoop stress on the hole boundary was the highest and it 
occurred in the bottom plate at the interface between the top and bottom 
plates at 90° to the load axis. The maximum stress concentration factor was 
found to be 5.7. Contact nonlinearity was confined to low levels of applied 
load (c„.< 30 MPa). At higher applied loads the response is almost linear. 
The radial stress distribution on the hole boundary can be represented by a 

<T r n 

cosine function of type = Cl COS 9 . Values of ‘a’ and ‘n’ depend 

on the through the thickness location on the hole boundary. 

(3) Elastic friction has negligible effect on local stresses (hoop) and hence it can 
be ignored. 

(4) For the clampup case, again the maximum stresses occurred on the hole 
boundary and in the bottom plate. The clampup force did decrease the 
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stresses at the critical location and its effect can be summarized in the form 
of a simple linear equation. For 10% clampup level, maximum hoop 
stress location changed from 75° at low load levels to 90° for high load 
levels. But for 25% and 35% clampup it remains at 75° for almost all load 
levels. The Clampup effect was confined to the region under the rivet head. 
Beyond the rivet head it’s effects were significantly reduced. 

(5) Interference is major factor that impacts the local stresses (hoop stress) 
around the rivet hole. Interference reduces the rate of increase of local 
stresses with remote loads. This causes the local (hoop) stresses to be 
lower than the neat fit results at high load levels. 

Elastic-Plastic Analysis of Two-Rivet Joint 

Elastic-plastic analysis of neat fit two-rivet lap joint was conducted. The material 
was assumed to be multi-linear stress-strain response. Von Mises yielding and isotropic 
hardening. The incremental elastic-plastic analysis was conducted to a remote stress of 
about 100 MPa. Hoop stresses at the critical hole at the joint was examined. The results 
showed that plastic deformation reduces stress concentration factor from 5.4 to 3.45, the 
hoop stress at 0 = 90° and z = 0 is nearly same as the maximum hoop stress and the 
bending stress reduce to nearly zero at large plastic deformation at the joint. 
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